High order Chin actions in path integral Monte Carlo 



K. Sakkos, J. CasuUeras, and J. Boronat 

Departament de Fisica i Enginyeria Nuclear, 
Universitat Politecnica de Catalunya, 
■ Campus Nord B^-BB, E-08034 Barcelona, Spain 

O ■ (Dated: March 16, 2009) 

^ ; Abstract 

5h ■ High order actions proposed by Chin have been used for the first time in path integral Monte 

Carlo simulations. Contrarily to the Takahashi-Imada action, which is accurate to fourth order only 
for the trace, the Chin action is fully fourth order, with the additional advantage that the leading 
I fourth and sixth order error coefficients are finely tunable. By optimizing two free parameters 

entering in the new action we show that the time step error dependence achieved is best fitted with 
' 5^ '■ a sixth order law. The computational effort per bead is increased but the total number of beads 

^ ■ is greatly reduced, and the efficiency improvement with respect to the primitive approximation is 

approximately a factor of ten. The Chin action is tested in a one-dimensional harmonic oscillator, 
a H2 drop, and bulk liquid ^He. In all cases a sixth-order law is obtained with values of the number 
^ I of beads that compare well with the pair action approximation in the stringent test of superfluid 

S '■ ^He. 

^ ■ PACS numbers: 31.15.Kb,02.70.Ss 
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I. INTRODUCTION 



Path integral Monte Carlo (PIMC) methods provide a fundamental approach in the 
study of interacting quantum many-body systems at low temperatures.^'^'- This accurate 
simulation tool relies on the well-known convolution property of the thermal density matrix 
which allows for estimating the density matrix at low temperature from their knowledge at 
higher temperature, where the system is well described by classical statistical mechanics.-i^ 
The partition function Z of the quantum system is then written as a multidimensional 
integral with a distribution law that resembles the one of a closed classical polymer with an 
inter-bead harmonic coupling. If one either ignores the quantum statistics of the particles 
(boltzmanons) , or they are bosons, then the statistical distribution law is positive definite and 
can be interpreted as a probability distribution function which can be sampled by standard 
Metropolis Monte Carlo methods. The mapping of this finite-temperature quantum system 
to a classic system of polymers was already suggested by Feynman^ and implemented in 
practice by Barker,- and Chandler and Wolynes^ in their pioneering works. 

In the simplest and most common implementations of the PIMC method the density 
matrix at high temperature, which constitutes the building block of the polymer chain, is 
considered fully classical, i.e., kinetic and potential parts of the action are splitted ignoring 
any contribution arising from the non-conmutativity of the kinetic K and potential V opera- 
tors. This approximation is known as primitive action (PA) and is accurate to order (/3/M)^, 
with (3 the inverse temperature and M the number of convolution terms (beads). When the 
temperature decreases the number of beads to reach convergence in terms of 1/M increases 
and the system becomes more and more quantum. In the classical limit, PIMC reduces to 
classic MC since the polymeric chain reduces to a single point (M = 1). On the other hand, 
near the zero-temperature limit the chains acquire extensions comparable to the size of the 
simulation box. PA is accurate to study semiclassical systems in which the quantum effects 
are comparatively small and therefore the number of beads in the asymptotic regime is low 
enough for an efficient sampling. However, if the interest lies in a fully quantum regime at 
very cold temperatures, M increases very fast making simulations hard, if not impossible, 
due to very low efficiency in the sampling of the long chains involved. 

The study of fully quantum fluids and solids require better actions than the PA. An 
accurate action to deal with very low temperatures, down to superfluid regimes, relies on the 
pair-product approximation^ in which the basic piece of the PIMC chain is the exact action 
for two isolated particles. The pair-product action has been used extensively in the study 
of superfluids and it is specially accurate for hard-spherelike systems such as liquid ^He.- 
However, its use for non-radial interactions is much more difficult due to the complexity 
of the pair density matrix. A different approach to improve the PA relies on the use of 
higher-order actions obtained directly from the exponential of the Hamiltonian, exp(— /3if). 
Working in this direction, Takahashi and Imada^'' and later on, and independently, Li and 
Broughtonii managed to write a new action (TIA) which is accurate to fourth order (/3/M)^ 
for the trace. To this end, they incorporated in the expression for the action the double 
commutator [[V^, K], V^]. More recently, Jang et al}"^ proposed a new action based on the 
Suzuki factorizationi^^ and accurate to fourth order. In spite of being a real fourth order 
action, and therefore better than the TIA, the results obtained by Jang et ai— do not show 
a significant improvement towards the reduction of the number of beads in the asymptote 
/3/M ^ respect to the TIA. 

In the present work, we introduce a new family of high order actions based on the sym- 
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plectic developments of Chini^ and Chin and Chen^^ that have already proved its efficiency 
in solving the Schrodinger equation,— problems in classical mechanics,-^^ and in the imple- 
mentation of evolution operators in density functional theory. ^'^ In Chin's factorization of the 
evolution operator all the coefficients are i) positive (and therefore directly implementable 
in Monte Carlo simulations) and ii) continuously tunable, which makes possible to force 
the error terms of fourth order to roughly cancel each other for the specific range of time- 
step values of interest for the actual simulations. We have introduced this new action in 
our PIMC algorithm and all the results presented in this article show that the accuracy 
achieved in terms of P/M is sixth order in practice in spite of the fact that the action is 
only accurate to fourth order. The attained reduction in the number of beads makes this 
algorithm competitive with the pair-product approximation in the study of quantum fluids 
and solids at ultralow temperatures. 

The rest of the work is organized as follows. The action and the method used for the 
PIMC simulation is contained in Sec. II. In Sec. Ill, we present results obtained with the 
Chin's action for the one-dimensional harmonic oscillator, a drop of H2 molecules, and bulk 
liquid ^He at low temperatures. Finally, a brief summary and the main conclusions are 
reported in Sec. IV. 



II. FORMALISM 

At finite temperature, the knowledge of the quantum partition function 

Z = Tre-^^ (1) 

allows for a full microscopic description of the properties of a given system, with (3 = 1/T 
and H = K + V the Hamiltonian. In a A^-particle system, the kinetic operator is given by 

i=l 

and the potential one, assuming pairwise interactions, by 

N 



V = J2v{r,,). (3) 



The non-conmutativity of the quantum operators K and V makes impractical a direct 
calculation of the partition function Z (Q. Instead, all practical implementations intended 
for Monte Carlo estimations of Z rely on its convolution approximation 

^-m+v) = fe-^iK+v)Y , (4) 



with e = P/M, where each one of the terms in the r.h.s. corresponds to a higher temperature 
MT. In the most simple approximation, known as primitive action (PA), the kinetic and 
potential contributions factorize 

e-^(i^'+v)^g~.A-g-.y ^ (5) 
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the convergence to the exact result being warranted by the Trotter formulap^ 



The primitive action is a particular case of a more general form in which one can decom- 
pose the action, that is 



e 



-£{K+V) 



Y^^-Uek ^-v^eV ^ (7) 



with {ti,Vi} parameters to be determined according to the required accuracy of the approx- 
imation. As we are interested in the Monte Carlo implementation of Eq. (j7]), all these 
parameters must be positive. However, the Sheng-Suzuki^^i^i theorem proves that this is 
impossible beyond second order in e. Recently, Chin^ has analyzed these expansions ([7]) for 
the estimation of the quantum partition function where only the trace is required ([T]). As 
the trace is invariant under a similarity transformation S, it could be possible to improve 
the order of the approximation by a proper choice of S. However, Chin^^ has proved that 
expansions like ([7]) can not be corrected with S beyond second order, generalizing in this 
form the Sheng-Suzuki theorem.— i2i 

In order to overcome the limitation imposed by the Sheng-Suzuki2Si2i theorem it is neces- 
sary to include in the operator expansion terms with double commutators such as [[V", i^], y]. 
In a recent work, Chin^^ has proved that a second-order algorithm (PA) can be corrected 
by a similarity transformation if that commutator is introduced. The result yields the TIA, 
with a trace that is accurate to fourth order. The TIA improves significantly the accuracy of 
the PA in PIMC simulations but not enough to deal properly with fully quantum fluids. In 
order to make a step further it is necessary to work directly with real fourth-order actions. 
Chin and Chen-^ have worked out a continuous family of gradient symplectic algorithms 
which are accurate to fourth order and that have proved to be extremely accurate in the 
resolution of classical and quantum problems. Now, we extend the focus of its applicability 
to the PIMC algorithm. 

The fourth-order action we have used is a two-parameter model given explicitely by^ 

hereafter referred as Chin action (CA). The generalized potentials W resemble the one that 
appears in the TIA since both incorporate the double commutator 



m 

i=l 



with Fi the force acting on particle i, 



N 



F, = Y. v.y{n,) . (10) 



The potentials W in Eq. ([8]) are explicitely 



Wa, = V+^a,e^(^f2\FA (11) 
m-2a. =V+^{1- 2a,)e' f- V \fA ■ (12) 
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The parameters in Eqs. ([8]) and fllllll2p are not all independent and can be written as a 
function of only two, ai and to, which are restricted to fulfill the conditions 



< ai < 1 (13) 
< < ^ ( 1 - 4^ ) • (14) 



2 V V3 

The rest of parameters are obtained from the two independent ones {ai,to} according to 
the equations 



1 - 2to 6(1- 2tof 



(15) 



Un = 

" 12 

Vl = — ^- (16) 

6(l-2to)2 ^ ' 

W2 = 1 - 2t;i (17) 
ti = \-t^. (18) 

The accuracy of the CA depends on the particular values of a\ and to that have to be 
numerically optimized. Each one modifies the action in different directions: to controls the 
weight of the different parts in which the kinetic part is splitted ([8]) and a\ the weight of 
each part in which the double commutator is divided (fT2l) . 

Restricting first our analysis to distinguishable particles, the quantum partition function 
Z ([1]) can be obtained through a multidimensional integral of the M terms (beads) in which 
it is decomposed, 

„ M 

Z= dRi... dRM n ' (19) 

a=l 

with R = {ri, . . . , tn} and Rm+i = Ri- In the rest of the work, Latin and Greek indexes 
are used for particles and beads, respectively. The density matrix of each step in Eq (fT9|) is 
written in the CA, 



p{Ra, -R( 



■a+l) — 

, 9N/2 / T \ 3Af/2 , 



0*2; I 9+2+ ; / dRaAdRaB exp ^ 
^ fl 1 1 \ 

N 

-e ^ {viV{ro,^,j) + V2V{r 

- e^uo — V (ai|F„,,|2 + (1 - 2al)|F,A,^|' + ai\F^B,i\^) } . (20) 
i=i J 

According to the CA, each elementary block of width e is splitted into three, with two middle 
points that we have denoted as A and B in the above expression for the density matrix fl20l) . 
Each one of the three internal beads resembles a bead in the TIA approximation in the sense 
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that the beads of the same type {a, aA, aB} interact through a generahzed Takahashi-Imada 
potential, but with different weights in front of the double-commutator term fllllll2p . The 
estimators for the total and partial energies of the system are therefore similar to the ones 
derived in the TIA approximation. 

The total and kinetic energy per particle can be readily derived from first derivatives of 
the quantum partition function Z, 



E 
N 
K 
N 



1 dZ 



NZ dp 
m dZ 



Nf3Z dm 



(21) 
(22) 



the potential energy being the difference V/N = E/N — K/N. The potential energy can 
also be derived through the relation^^ 



0(R) 



dZiV + AO) 



PZ{V) 



(23) 



A=0 



that can be used also for the estimation of other coordinate operators. From the definition 
( I22ll . which is known as thermodynamic estimator, the kinetic energy per particle results in 



K 



th 



1 



m 



N 
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2e MN \2h^e'^ 



Tlrj,-—e\oW, 



m 



MN 



(24) 



with 



rpt 



M N 

EE 

a=l 1=1 



faA,i 



+ — [raA,^ 
tl 



raB,i 



+ — [r^B, 



'a+l,i) 



(25) 



and 



M N 

W^M7v = 5ZZlKl^--il' + (l-2ai)|F,A,ir + ai|F„B,,|2] , (26) 

a=l i=l 

The potential energy can be calculated from the difference between the total energy and 
or by means of the general relation (125]) with identical result. 



the kinetic one 



V_ _ _J_ 

N ~ MiV 



Vj 



MN 



2 — e^oWmn 
m 



(27) 



with 



M N 

Vmn = X] + ^2V^(w,y) + viV{raB,ij)] , (28) 

a=l i<j 

and Wmn given by Eq. (1261) . This term that appears both in the total and kinetic energy, 
but with different weight, comes from the derivative with respect to f3 and with respect to 
m, respectively, of the W potentials fllllll2p . which depend explicitely on the temperature 
and the mass. 

The variance of the thermodynamic estimation of the kinetic energy fl2^ can be rather 
large and increases with the number of beads M.—^ This well-known problem is generally 
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solved by using the virial estimator^- of the kinetic energy which rehes on the invariance 
of the partition function under a scahng of the coordinate variables r Xr. The centroid 
version of the virial estimator for the CA is given by 



K""" 3 1 /I n'^ 



TT:^ + 7777 o ^Miv + - e'uoiWMN + Ymn) • (29) 



N 2(3 MN \2 ''''' m 
In this expression, Wmn is given by Eq. (I26l) and the two other terms are explicitely 



M N 



Tmn = X] X] [Vlir a,i - To^i) F^^i + V2{raA,i - ^'o.i) + Vi{raB,i - ro,i) FaB,i\ (30) 

a=l i=l 

and 

M N N 
a=l i=l j^i 

+ (1 - 2al){r^A,^ - ro,iTT{aA,^,j)''^{F^A,^ - F^Aj)b 

+al{r^B,^ - ro,iTTiaB,i,jf^{F^B,i - F^B,j)b] • (31) 

In the above expressions, r,, j is the center of mass (centroid) of the chain representing the 
atom i. The indexes a, b in the definition of Y^n dSI]) stand for the Cartesian coordinates 
and an implicit summation over repeated indices is assumed. The tensor T appears also in 
the virial estimation of the kinetic energy in the TIA approximation and is explicitely given 
by^ 



■ ■ dv^ ■ ■ 



(32) 



where 7 stands for the three different types of internal beads, a, aA^ and a-B, and 5^ = 1 if 
a = h and otherwise. 

The implementation of the CA in the PIMC algorithm is similar to the one for the 
TIA. Going from TIA to CA, one has to split a single bead into three (with different link 
lengths), but in these new beads the different atoms interact with a similar generalized 
potential, coming from the double commutator. A simple inspection on the equations of 
the density matrix and energies for the CA and TIA shows that the complexity of both 
actions are essentially the same and require, for example, the same order of derivatives of 
the interatomic potential V{t\ 

A final concern that any PIMC calculation must afford is whether the sampling method 
has the necessary efficiency in the movement of the chains to avoid the slowing down that 
can appear for long chains when using only individual bead movements. In this sense, the 
implementation in the algorithm of collective smart movements is crucial. To this end, we 
use the staging method^i^i2^ combined with movements of the center of mass of each one 
of the atoms. In the CA, the length of each chain is not the same and therefore one has to 
generalize the staging method. This generalization is discussed in Appendix A. 
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III. RESULTS 



We have studied the accuracy of the high-order action proposed by Chin in three different 
systems: a one-dimensional harmonic oscillator, a drop of H2 molecules, and bulk liquid ^He. 
The harmonic oscillator is a very simple system but it has the advantage of its exact analytic 
solution at any temperature which allows for an accurate test of the method. The H2 drop 
is a more exigent test and it has been adopted in previous works as benchmark for testing 
different actions. Finally, liquid ^He is also a well-known system with a wide experimental 
information available for comparison and with the appealing feature of remaining liquid even 
at zero temperature. In all the simulations we have used a sampling method which combines 
movements of the center of mass of the chains and smart movements of chain pieces with 
the staging technique .->2^>2& The length of the staging chain and the maximum value of the 
displacement of the center of mass are chosen to achieve an acceptance rate of 30%-50%. 
The results presented below have been checked to be stable with respect to the frequency 
and size of both movements. 



A. Harmonic oscillator 

In our first application, we consider a particle in a one-dimensional harmonic oscillator 
with Hamiltonian 

H = — + -mu^x^ . (33) 

As it is well know, this problem can be exactly solved at any temperature and, in particular, 
the energy is given by^ 

E = ^hujcoth{(3hu/2) . (34) 

In the PIMC simulations we have taken uj = h = m = 1 and the results presented correspond 
to temperatures T = 0.1 and 0.2. The energies obtained at T = 0.2 with different actions 
and as a function of the number of beads M are contained in Table [H and have to be 
compared with the exact value 0.50678. In order to reproduce the exact energy with these 
five digits PA requires the use of M = 512 and TIA of a quite smaller number M = 32.— 
This asymptotic value is reached in the CA case with a sizeable smaller number: M = 6 
and M = 5 for ai = and ai = 0.33, respectively. 

In Fig. [H, the energies for the different actions at T = 0.2 are plotted as a function of e. 
The lines on top of the PIMC data correspond to polynomial fits of the form 

E = Eo + Ase^ (35) 

to the energies close to the common asymptotic value £"0. From previous work— it is known 
that 6 = 2 and 4 for PA and TIA, respectively. In the case of the CA approximation 
the departure from £"0 is effectively of sixth order 5 = 6 in spite of the fact that CA is 
rigorously a fourth-order action. This result can be understood as a partial but effective 
cancellation between the leading errors of fourth and higher orders due to the fact that 
the CA ([8]) is written in terms of two parameters to and oi that can be optimized within 
some constraints fll3lll4p . We show in Fig. [2] the characteristic behavior of the energy as a 
function of e and for different choices of to and a fixed value for Oi. The lines are fits to the 
PIMC data for values of to in the range 0.9-0.15 (top to bottom for large e values in the 
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M 




-E-TIA 


EcA{ai = 0) 


Eca(«i = 0.33) 


2 


0.30755 


0.44702 


0.50444 


0.50643 


3 






0.50649 


0.50675 


4 


0.43162 


0.50053 


0.50673 


0.50677 


5 






0.50677 


0.50678 


6 






0.50678 


0.50678 


7 

8 


0.48424 


0.50630 


0.50678 




16 


0.50085 


0.50675 






62. 


0.50528 


0.50678 






64 


0.50641 


0.50678 






128 


0.50669 








256 


0.50676 








512 


0.50678 









TABLE I: PA (E'pa), TIA (£^tia), and CA {Eqa) results for the one-dimensional harmonic oscil- 
lator at T = 0.2. 

0.5069 



0.5067 



W 0.5065 



0.5063 



0.5061 

0.0 0.5 1.0 1.5 2.0 2.5 

8 

FIG. 1: PIMC energy of a particle in a one-dimensional harmonic well as a function of e. Triangles, 
diamonds, squares and circles stand for PA, TIA, CA (ai = 0),and CA(ai = 0.33), respectively. 
The lines correspond to polynomial fits (j35p to the data. The errorbars are smaller than the size 
of the symbols. 

figure). When < 0.13 the asymptotic exact value is approached from above and contrarily 
from below when to ^ 0.13. By adjusting in a proper way the value of to it is therefore 
possible to achieve a nearly flat dependence with e and consequently to improve empirically 
the fourth-order accuracy of the CA. 

The estimation of the optimal value of the parameter to for a fixed ai can be better made 
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0.50680 



0.50679 



W 0.50678 



0.50677 - 



0.50676 




0.0 0.2 0.4 0.6 0.8 



FIG. 2: Departure from the asymptotic energy £"0 for different values of to and a fixed value for ai 
(01 = 0.33). Prom top to bottom, to = 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, and 0.15. The temperature 
is T = 0.2. 




0.48 



0.04 0.08 0.12 0.16 0.2 

to 



FIG. 3: Energies as a function of the parameter to for ai = 0.33, and for different number of beads, 
at T = 0.1. Squares, M = 2; filled squares, M = 3; circles, M = 4; filled circles, M = 5; up 
triangles, M = 6; up filled triangles, M = 7; down triangles, M = 8. 
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0.50006 




FIG. 4: Energies as a function of e for different values of the parameter ai at T = 0.1. Squares, 
ai = 0; filled squares, oi = 0.14; circles, ai = 0.25; filled circles, ai = 0.33; triangles, ai = 0.45. 
The optimal values of to are, respectively, 0.1430, 0.0724, 0.1094, 0.1215, and 0.1298. The lines are 
polynomial fits (j35p to the PIMC data. 



by computing the energies for an increasing number of beads, starting for a low value. The 
behavior observed, which is similar for any allowed ai value, it is shown in Fig. [3] for the 
particular case of ai = 0.33 and T = 0.1. The lines on top of the data are guides to the 
eye and each one of them correspond to a particular M in the range M = 2-8. As one can 
see, there are two small intervals in the to scale where the curves tend to intersect (around 
0.12 and 0.16) and therefore where the convergence to the asymptotic value is faster. By 
examining the energies inside these two regimes one observes that the dependence with M 
is slightly smoother in the first one corresponding to smaller to's- Thus, we work in this 
first regime and normally by selecting a value in which the approach to the exact energy 
is from below. For example, for this value ai = 0.33 we have chosen to = 0.1215. It is 
worth noticing that this optimal value does not depend on temperature and consequently 
it can be adjusted working at higher temperatures where the number of beads required to 
achieve convergence is always smaller and thus cheaper from a computational point of view. 
We have verified that the best value of to, for a given ai, obtained through this numerical 
optimization agrees with the analytical relation for the harmonic oscillator obtained in Ref. 



271 , which predicts the optimal parameters that cancel exactly the fourth-order error terms. 

The optimal to depends on the particular value of Oi and the achieved accuracy depends 
also slightly on ai. This is explicitly shown in Fig. H] where the dependence of the energy 
on 6 is plotted for values of ai ranging from to 0.45. The results correspond to T = 0.1 
and the lines correspond to fits fl35l) with 6 = 6. In all cases the accuracy is of the same 
order but the best performance is achieved for ai = 0.33 which, according to the expression 
of the CA, is the case where the generalized potential W f fTT1[T2l) acts with the same weight 
in the three points in which a single step e is splitted. 



11 



M 


(i?/iV)pA 


(i?/iV)TIA 




2 






-40.44(5) 


4 






-28.77(3) 


8 


-45.28(3) 


-31.17(3) 


-21.27(2) 


10 








12 






-19.13(3) 


14 








16 


-30.59(3) 


-23.48(3) 


-18.32(2) 


18 








20 






-17.95(3) 


24 






-17.81(2) 


28 






-17.73(2) 


32 


-22.97(3) 


-19.49(3) 


-17.66(2) 


36 






-17.68(2) 


48 






-17.68(2] 


64 


-19.57(3) 


-18.05(3) 




128 


-18.28(3) 


-17.78(3) 




256 


-17.89(3) 


-17.72(3) 




512 


-17.76(3) 







TABLE II: PA, TIA, and CA results for the energy per particle of a drop composed by = 22 H2 
molecules at T = 6 K for different values of M. All the energies are in K. The figures in parenthesis 
are the statistical errors affecting the last digit. 

B. H2 drop 

The case study of a drop composed by a few number of hydrogen molecules has been 
used in the past to compare the efficiencies of different PIMC methods and, in particular, 
of several approximations for the action. For this purpose, this system was used for the 
first time by Chakravarty et al?^ to compare the efficiency of Fourier vs. standard PIMC 
methods. Later on, it was studied by Predescu et alM in a comparative analysis of energy 
estimators and by Yamamoto^ in a fourth-order calculation of small atomic and molecular 
drops. 

The drop studied is composed by = 22 H2 molecules which are considered spherical 
since we restrict our calculation to the J = state, i.e., to parahydrogen. The interaction 
potential is of the form 

N N 

V{n, . . . , r^) = ^ V^ir,,) + J2 K(r.) , (36) 

i<j i=l 

with V2 the intermolecular interaction, assumed to be of Lennard- Jones type 



cr 



12 



-18 



^ -20 - 



-22 



-24 








0.004 



0.008 

e(K"') 



0.012 



0.016 



FIG. 5: Energy of the H2 drop at T = 6 K as a function of e for different actions: PA, diamonds; 
TIA, squares; circles, CA. The parameters of the CA are ai = and to = 0.175. The hnes are 
polynomial fits (|35p to the PIMC data. 



and K is a confining potential introduced to suppress possible evaporation of molecules, 

20 



CM 



Re 



(38) 



In Vci Rem = X]j=i ^i/^ is the position of the center of mass of the drop and Rc controls 
the maximum allowed distance of a particle to -Rcm- As in previous work in this problem,— 
we have chosen R^ = 4cr, with Lennard- Jones parameters a = 2.96 A and e = 34.2 K fl37|) . 

In Table HIl we report our results for the energy of the drop using different number of 
beads M and several approximations for the action: PA, TIA, and CA. The asymptotic 
(unbiased) energy is obtained in the limit e — > (M 00): this is achieved with M ~ 512, 
128, and 32 for PA, TIA, and CA, respectively. Therefore, also in this more exigent test the 
CA shows its appreciably higher efficiency with respect to TIA and other published data 
of the same problem obtained with Suzuki actions.— The value of e required to reach the 
asymptote in CA is comparable with the one observed in a previous study of this H2 drop 
using the pair action approximation.— As it is shown in Fig. [5l the e — behavior is well 
reproduced by the polynomial fit (!35l) with different exponents: 6 = 2, 4, and 6 for PA, 
TIA, and CA, respectively. Our best result for the energy of the drop is obtained from the 
simulation with the CA: the total energy is E/N = —17.68(2) K, and the potential and 
kinetic energies are V/N = —47.82(3) K and K/N = 30.14(2) K , respe ctively. These three 
energies are in agreement with the previous estimations of Refs. 28 29 30l. 

Going down in temperature, the PIMC calculation becomes harder due to the increase 
in the number of beads necessary to achieve convergence and to the simultaneous decrease 
of the acceptance rate in the staging movements. In order to be effectively able to compute 
properties in these ultracold regimes the action in use has to be accurate enough to reduce 
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£(K"') 

FIG. 6: Energy of the H2 drop at T = 1 K as a function of e using the CA. The hne is a polynomial 
fit ([35]) to the PIMC data with 5 = Q. 

M to a manageable level. To this end, we have computed the properties of the N = 22 drop 
at T = 1 K. In this simulation, as in the one at T = 6 K, we have ignored the Bose statistics 
of the molecules. The dependence of the total energy with e is shown in Fig. O As one can 
see, the value of e (M = 180) in the asymptote is approximately a factor of two smaller than 
the one at T = 6 K (Fig. [5]) but a sixth-order behavior is again observed (solid line in Fig. 
[6]). The energies in the asymptotic regime are E/N = —22.47(5) K, V/N = —52.7(2) K,and 
K/N = 30.2(2) K. It is worth noticing that at this lower temperature PA would require the 
use of M ~ 3000 and TIA of M ~ 1200, values which make quite unreliable their use in the 
deep quantum regime. 

C. Liquid ^He 

As in our previous work,24 we have studied the accuracy of the CA in a fully many-body 
calculation, deep in the quantum regime, as it is liquid ^He. We consider a bulk system at 
a density p = 0.02186 A"'^ and at two temperatures, T = 5.1 and 0.8 K. The calculation is 
performed within a simulation box of 64 atoms with periodic boundary conditions and with 
an accurate Aziz potential. '^^ In order to correct for the finite size of the system we have 
added standard potential energy tail corrections relying on the assumption that beyond L/2, 
with L the size of the box, the medium is homogeneous, i.e. g{r) ~ 1. 

Table IIIII contains PIMC results for the energy at T = 5.1 K using different number 
of beads and for several actions: PA, TIA, and CA with ai = and ai = 0.33. The 
parameter to has been optimized for both values of oi: to = 0.17 and 0.082 for ai = 
and 0.33, respectively. As we noted in our previous study,— the introduction of the double 
commutator within the TIA reduces sizably the number of beads with respect to PA: the 
value M = 256 of PA turns to M = 64 for the TIA. The latter M is again considerably 
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M 


(i?/iV)pA 


(i?/iV)TiA 


{E/N)cK (ai = 0) 


{E/N)cA («i = 0.33) 


3 






-5.59(3) 


-5.48(3) 


4 






-4.51(3) 


-4.40(3) 


5 






-3.91(3) 


-3.77(3) 


6 






-3.53(3) 


-3.43(3) 


8 


-9.29(3) 


-6.30(3) 


-3.16(3) 


-3.08(3) 


10 






-3.01(3) 


-2.89(3) 


12 






-2.91(3) 


-2.83(3) 


14 






-2.86(3) 


-2.81(3) 


20 






-2.81(3) 




16 


-5.85(3) 


-3.97(3) 






32 


-3.98(3) 


-3.04(3) 






44 




-2.93(4) 








-o.iy(4j 


-z.oo(4j 






100 




-2.83(4) 






128 


-2.92(4) 


-2.82(4) 






256 


-2.84(4) 








512 


-2.81(4) 









TABLE III: PA, TIA, and CA (oi = and oi = 0.33) results for the energy per particle of bulk 
liquid ^He at T = 5.1 K for different values of M. All the energies are in K. The figures in 
parenthesis are the statistical errors affecting the last digit. 



reduced by using the CA since the convergence to the exact energy is reached for a value as 
low as M = 12 (oi = 0.33). The different dependence on e of the three models for the action 
is shown in Fig. [7|at T = 5.1 K. We can observe that the departure from the asymptote Eq 
follows the same power-law dependence ( l35i) than in the harmonic oscillator and H2 drop 
previously analyzed: 5 = 2, 4, and 6 for PA, TIA, and CA, respectively. 

The accuracy of the CA has been a bit more stressed by repeating the PIMC simulation 
at a lower temperature, T = 0.8 K. At this temperature, the quantum effects are bigger 
than at T = 5.1 K and the use of TIA, and even more PA, is completely unreliable due to 
the large number of beads that are necessary to eliminate the bias due to a finite M value. 
At 0.8 K, liquid "^He is below the A transition (Tx = 2.17 K) and therefore it is superfluid 
making absolutely necessary the sampling of permutations to accomplish with its boson 
statistics. The exchange frequency is drastically reduced over T\ and thus the inclusion 
of permutations at the higher temperature (5.1 K) does not modify the results presented 
above. Results for the energy per particle at T = 0.8 K as a function of e and using the 
CA are shown in Fig. [HI The simulations have been carried out including or not the right 
symmetry of the thermal density matrix; as one can see, and is well known, the inclusion of 
permutations leads to a decrease of the energy which is in agreement with the experimental 
data on this system. The sampling of permutations is performed within the widely used 
method proposed by Pollock and Ceperley.-'^ The accuracy of the CA is the same including 
or not permutations in the sampling , i.e., sixth order in e (solid lines in Fig. [S]). The values 
of e at which the asymptotic trend is observed are similar to the ones achieved using the 
pair action,^ which is the most widely used approximation to deal with superfiuids within 
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FIG. 7: Energy per particle of liquid ^He at T = 5.1 K and density p = 0.02186 as a function 
of e using several actions: PA, squares; TIA, diamonds; CA (ai = 0), circles; CA (ai = 0.33), 
triangles. The lines are polynomial fits (j35p to the PIMC data with (5 = 2, 4, and 6 for PA, TIA, 
and CA, respectively. 
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FIG. 8: Energy per particle of hquid ^He at T = 0.8 K and density p = 0.02186 A'^ function 
of e using the CA. Circles and squares correspond to a simulation including symmetrization of the 
density matrix or not, respectively. The lines are polynomial fits (|35p with 6 = 6. 
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CPU cost per bead 


Decrease of M 


Efficiency 


PA 


1.0 


1 


1.0 


TIA 


2.9 


4 


1.4 


CA 


7.2 


58 


8.0 



TABLE IV: Comparison among the efficiencies of the PA, TIA, and CA in PIMC. 



the PIMC formalism. 

IV. CONCLUSIONS 

In the last decades there has been a continued effort for improving the action to be used in 
PIMC simulations beyond the PA. Working directly on the exponential of the Hamiltonian, 
Takahashi and Imada^'' introduced in the action the double commutator K], V^] and 
showed that the new algorithm (TIA) was of fourth order. As showed in previous work,— 
the TIA reduces significantly the number of beads to reach the asymptotic limit and therefore 
it can be very useful in quantum systems if the temperature is not very small. However, if 
one is interested on achieving lower temperatures, deep in the quantum regime, the TIA is 
still not accurate enough since the number of beads required is yet too large. As pointed 
out by Chin,— this relative failure of TIA is due to the fact that the TIA is in fact a fourth- 
order action but only for the trace. Posterior attempts of improving the action,— based 
on actions proposed by Suzuki, did not show a significant enhancement of the efficiency of 
the method with respect to TIA. 

Following with the aim of going a step further on the improvement of the action for PIMC 
applications, we have used for the first time the new developments of Chin^^ii^ that have led 
to full fourth-order expansions of the exponential of the Hamiltonian. The resulting action 
(CA) is more involved than the TIA and Suzuki action but still it is readily implementable 
starting on a TIA approach since the basic ingredients of the CA are already contained in 
the TIA. In Table llVt we compare the efficiency of the three models for the action; the 
numbers correspond to the calculation of liquid "^He but are similar to the ones obtained 
in the study of the H2 drop. Considering as 1 the CPU time per bead and the efficiency 
in the PA one observes that the cost per bead in TIA is 2.9 and the number of beads is 4 
times smaller, leading to en efficiency 4/2.9 = 1.4. In the CA, the CPU cost per bead is 
larger since every step e is splitted into three but the big decrease in the number of beads 
(58 times smaller than in PA) results in an efficiency 8 times larger than the PA. 

In the Chin action there are two parameters (to, cti) that have to be optimized. The 
search of these parameters for a given system is not difficult since they are independent of 
temperature and therefore they can be found at higher temperatures where the number of 
beads is very small and consequently the simulations very short in CPU time. As we have 
discussed in the harmonic oscillator problem, the exact energy is crossed by changing Iq 
because the departure from the asymptote can be upwards or downwards depending on its 
particular value. This behavior makes possible to search for optimal parameters that improve 
the efficiency of the action from fourth-order to an effective sixth order. This behavior is 
not only characteristic of the harmonic oscillator but completely general as we have verified 
in real many-body problems as the H2 drop and bulk liquid ^He. This sixth-order law 
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is maintained even at superfluid temperatures making that the number of beads required 
for the simulations is completely manageable. Therefore, the CA is a realistic alternative 
to the pair action for dealing with quantum liquids and solids in the superfluid regime at 
temperatures close to zero. Its implementation is not much more involved than the TIA, 
easier to use than the pair action, and useful also for problems with non-radial interactions 
where the application of pair action is much more involved. Work is also in progress to 
use the CA in a path integral ground state (PIGS) approach to study the limit of zero 
temperature, the initial simulations showing also a sixth-order accuracy and convergence to 
the exact energies with a few number of beads. 
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APPENDIX: STAGING TRANSFORMATION 



The staging technique allows for a direct sampling of the free (kinetic) part of the action, 
i.e.. Metropolis test is not necessary. In this Appendix, we generalize the standard staging 
method to the one required for the CA action where the width of a given bead depend of its 
type {tie, 2toe) fl20|) . If two points of the chain representing an atom are considered fixed, 
Vq and Vm, one is interested in transforming the free action between these two extremes 



S = exp 



M 



-E 



(A.l) 



into the staging one 



M-l 



'S'st = C(ro, tm) exp <^ - ^ g„[r„ - (a„rc,_i + h^rM)]^ > • (A. 2) 



a=l 



The constant C depends only on the fixed positions Vq and tm and therefore it is not 
important for the purpose of the sampling. By imposing that both actions (1A.1IIA.2P are 
equal one can derive sequential relations for the staging coefficients g^, Cq,, and ha as a 
function of the known ones of the original action Cq, 

la = + — (la+i'^'ii+i (^-3) 

aa = Ca/qa (A.4) 
ba = {Qa+lO-a+lba+l + CMSa,M-l) /la , (^-5) 

with Sa,M-i = lifa = M — 1 and otherwise. These relations have to be applied recursively 
from M — 1 to 1 with the starting conditions Qm = clm = &m = 0. 
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Once all the staging coefficients are determined through Eqs IA.3tlA.5l the new positions 
Va are recursively obtained from 1 to M — 1 using Gaussian displacements. 
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